── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ lubridate 1.9.3 ✔ tibble 3.2.1
✔ purrr 1.0.2 ✔ tidyr 1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ purrr::flatten() masks jsonlite::flatten()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dbscan)
Attaching package: 'dbscan'
The following object is masked from 'package:stats':
as.dendrogram
library(mclust)
Package 'mclust' version 6.1.1
Type 'citation("mclust")' for citing this R package in publications.
Attaching package: 'mclust'
The following object is masked from 'package:purrr':
map
library(randomForest)
randomForest 4.7-1.2
Type rfNews() to see new features/changes/bug fixes.
Attaching package: 'randomForest'
The following object is masked from 'package:dplyr':
combine
The following object is masked from 'package:ggplot2':
margin
The function kmeans() performs \(K\)-means clustering in R. We begin with a simple simulated example in which there truly are two clusters in the data: the first 25 observations have a mean shift relative to the next 25 observations.
The \(K\)-means clustering perfectly separated the observations into two clusters even though we did not supply any group information to kmeans(). We can plot the data, with each observation colored according to its cluster assignment.
#par(mfrow = c(1, 2))plot(x, col = (km.out$cluster +1),main ="K-Means Clustering Results with K = 2",xlab ="", ylab ="", pch =20, cex =2)
Here the observations can be easily plotted because they are two-dimensional. If there were more than two variables then we could instead perform PCA and plot the first two principal components score vectors.
In this example, we knew that there really were two clusters because we generated the data. However, for real data, in general we do not know the true number of clusters. We could instead have performed \(K\)-means clustering on this example with \(K=3\).
plot(x, col = (km.out$cluster +1),main ="K-Means Clustering Results with K = 3",xlab ="", ylab ="", pch =20, cex =2)
When \(K=3\), \(K\)-means clustering splits up the two clusters.
To run the kmeans() function in R with multiple initial cluster assignments, we use the nstart argument. If a value of nstart greater than one is used, then \(K\)-means clustering will be performed using multiple random assignments in Step~1 of Algorithm 12.2, and the kmeans() function will report only the best results. Here we compare using nstart = 1 to nstart = 20.
Note that km.out$tot.withinss is the total within-cluster sum of squares, which we seek to minimize by performing \(K\)-means clustering (Equation 12.17). The individual within-cluster sum-of-squares are contained in the vector km.out$withinss.
We strongly recommend always running \(K\)-means clustering with a large value of nstart, such as 20 or 50, since otherwise an undesirable local optimum may be obtained.
When performing \(K\)-means clustering, in addition to using multiple initial cluster assignments, it is also important to set a random seed using the set.seed() function. This way, the initial cluster assignments in Step~1 can be replicated, and the \(K\)-means output will be fully reproducible.
Let’s try k-means clustering on a more complicated dataset. Open questions:
How many clusters are in this dataset? (Hyperparameter choice)
How many can we expect k-means to detect?
df_hier %>%ggplot(aes(x = xpo, y = ypo)) +geom_point()
We can attempt to tune our clustering algorithm by observing the total within-cluster sum of squares over each level of the hyperparameter k. Keep in mind this is by no means a “correct” partition of the dataspace but rather a tool to aid in decisionmaking.
tune_kmeans =function(data, kmax) { result =c()for (i inc(1:kmax)) { km.out <-kmeans(data, i, nstart =20) result[i] = km.out$tot.withinss } df_viz =data.frame(k =c(1:kmax),within_ss = result)}df_viz =tune_kmeans(df_hier, 10)df_viz %>%ggplot(aes(x = k, y = within_ss)) +geom_line() +scale_x_continuous(breaks=seq(1,10,1))
The hclust() function implements hierarchical clustering in R. In the following example we use the data from the previous lab to plot the hierarchical clustering dendrogram using complete, single, and average linkage clustering, with Euclidean distance as the dissimilarity measure. We begin by clustering observations using complete linkage. The dist() function is used to compute the \(50 \times 50\) inter-observation Euclidean distance matrix.
The second argument to cutree() is the number of clusters we wish to obtain. For this data, complete and average linkage generally separate the observations into their correct groups. However, single linkage identifies one point as belonging to its own cluster. A more sensible answer is obtained when four clusters are selected, although there are still two singletons.
plot(hclust(dist(df_hier), method ="complete"),labels =FALSE, hang =-1, ann =FALSE, xlim =c(0,10))
set.seed(4)hc.complete <-hclust(dist(df_hier), method ="complete")df_viz = df_hierdf_viz$cluster =factor(cutree(hc.complete, 9))df_viz %>%ggplot(aes( x = xpo, y = ypo, color = cluster)) +geom_point()
Using hierarchical clustering we have a nice understanding of the hierarchical nature of our data. The dendrogram shows clearly the order in which the communities are constructed and can be used to detect nested communities. We can observe that the model is quite susceptible to the random noise between the “expected” clusters.
DBSCAN
If we use a density based approach the result changes slightly: The model fails in helping us understand the hierarchical data structure on the right. Keep in mind that the algorithm directly uses the density of observations in the dataspace. A mix of high and low clustering makes a proper choice for epsilon difficult.
It is however much less susceptible to noise as the minPts parameter does not allow the model fit arbitrarily small clusters, instead all the noise is grouped into one additional cluster.
eps <-5minPts <-3db <-dbscan(df_hier, eps = eps, minPts = minPts)df_viz = df_hierdf_viz$col =factor(db$cluster)df_viz %>%ggplot(aes( x = xpo, y = ypo, color = col)) +geom_point()
Model Based Clustering - Gaussian Mixture Model
Gaussian-mixture model is as the name implies a model-based approach, meaning we can get uncertainty estimations for each data point and each cluster!
The “classification” vector stores the highest probability class of each observation. The “uncertainty” vector stores the related uncertainty of that classification.
gmm =Mclust(df_hier)df_viz = df_hierdf_viz$col =factor(gmm$classification)df_viz %>%ggplot(aes( x = xpo, y = ypo, color = col)) +geom_point()
Unfortunately we loose some flexibility with a model based approach: The gaussian assumption does not hold for many real world datasets and we cannot really detect arbitrarily shaped communities as the model expects convex / ellipsoid distributions.
gmm =Mclust(df_snake)df_viz = df_snakedf_viz$col =factor(gmm$classification)df_viz %>%ggplot(aes( x = xpo, y = ypo, color = col)) +geom_point()
Compare this to the previous models: Hierarchical clustering identifies our expected communities when using single linkage and DBSCAN can be tuned to do so as well.
set.seed(4)hc.complete <-hclust(dist(df_snake), method ="single")df_viz = df_snakedf_viz$cluster =factor(cutree(hc.complete, 4))df_viz %>%ggplot(aes( x = xpo, y = ypo, color = cluster)) +geom_point()
eps <-5minPts <-3db <-dbscan(df_snake, eps = eps, minPts = minPts)df_viz = df_snakedf_viz$col =factor(db$cluster)df_viz %>%ggplot(aes( x = xpo, y = ypo, color = col)) +geom_point()
Application - Young People Survey:
Let’s use our newfound skills to classify some people. Below we find survey data of 1000 people in Slovakia aged 15 to 30. We will try to detect distinct communities using clustering.
First we load the dataset and keep only numeric data (we could also transform the categorical)
set.seed(4)hc.complete <-hclust(dist(df_youth), method ="complete")plot(hclust(dist(df_youth), method ="complete"),labels =FALSE, hang =-1, ann =FALSE)
df_viz = df_youthdf_viz$cluster =factor(cutree(hc.complete, 4))df_viz %>%ggplot(aes( x = Age,y = Number.of.friends, color = cluster)) +geom_jitter()
This tells us very little as we have no idea which variables to visualize! There appear to be at least two interesting communities here but to learn more we need to go back to supervised learning.
Warning in randomForest.default(m, y, ...): The response has five or fewer
unique values. Are you sure you want to do regression?
Interestingly enough we see factors stereotypically associated with “outsiders” as most relevant for belonging in the second cluster. We will continue here next week with PCA.
